Matter- wave bright solitons in spin-orbit coupled Bose-Einstein condensates 
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We study matter-wave bright solitons in spin-orbit (SO) coupled Bose-Einstein condensates (BECs) with 
attractive interatomic interactions. We use a multiscale expansion method to identify solution families for 
chemical potentials in the semi-infinite gap of the linear energy spectrum. Depending on the linear and spin-orbit 
coupling strengths, the solitons may resemble either standard bright nonlinear Schrodinger solitons or exhibit a 
modulated density profile, reminiscent of the stripe phase of SO-coupled repulsive BECs. Our numerical results 
are in excellent agreement with our analytical findings, and demonstrate the potential robustness of such solitons 
for experimentally relevant conditions through stability analysis and direct numerical simulations. 
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Introduction. Gauge fields are ubiquitous in physics, as 
they are relevant to the interactions of charged particles with 
electromagnetic fields III] or to fundamental interactions in el- 
ementary particle physics [2] . Ultracold atomic gases are con- 
sidered as an excellent candidate where a variety of artificial 
gauge fields can be realized; see Ref. [3] for a review. Such 
gauge fields have been recently studied in experiments 0, 
with binary Bose-Einstein condensates (BECs). Importantly, 
synthetic magnetic fields can produce spin-orbit (SO) interac- 
tions in a BEC consisting of (predominantly) two hyperfine 
states of 87 Rb, coupled by a Raman laser |5|]. 

SO-coupled BECs with repulsive interactions have become 
a topic of intense investigations. Different studies have re- 
vealed the existence of a "stripe phase" (consisting of a lin- 
ear combination of plane waves) (fl] and phase transitions be- 
tween it and states with a single plane wave or with zero mo- 
mentum IlZD - The existence of topological structures, such as 
vortices with ^\ or without rotation, Skyrmions [10] and 
Dirac monopoles full , as well as self-trapped states (solitons) 
of an effective nonlinear Dirac equation (NLDE), was also il- 
lustrated 11211 . While the above studies refer to BECs with re- 
pulsive interactions, to the best of our knowledge, SO-coupled 
BECs with attractive interactions have not been studied so far. 
The latter, is the theme of the present work. 

As it is known, attractive BECs can become themselves 
matter- wave bright solitons [13], i.e., self-trapped and highly 
localized mesoscopic quantum systems that can find a vari- 
ety of applications [14]. Here, we demonstrate the existence, 
stability and dynamics of matter-wave bright solitons in SO- 
coupled attractive BECs. In particular, starting from the corre- 
sponding mean-field model, we consider the nonlinear waves 
emerging in the semi-infinite gap of the linear spectrum. Sim- 
ilarly to the repulsive interaction case of Ref. |7|], we find three 
distinct states having: (a) zero momentum, (b) finite momen- 
tum, +fc or — ko, and (c) stripe densities formed by the inter- 
ference of the modes with ±fco momentum. We analytically 
identify these branches, in very good agreement with our nu- 
merical computations, and determine their spin polarizations. 
We also analyze the stability of these solutions, illustrating 
that branches (a) and (c) are generically stable, while branch 
(b) is stable for sufficiently small atom numbers. Hence, these 
newly emerging matter-wave solitons in SO-coupled BECs 



may be well within experimental reach. 

Model We consider SO-coupled BECs confined in a quasi- 
1D parabolic trap, with longitudinal and transverse frequen- 
cies uj x <C cj_l- In this setting, and for equal contributions of 
Rashba [15] and Dresselhaus [16] SO couplings (as in the ex- 
periment of Ref. [5]), the mean-field energy functional of the 
system is E = £dx, with: 

e=±{&H 9 + 9nW 4 + 922\i>i\ 4 + 2 ffl2 |V t | 2 |V';| 2 ),(l) 

where \l> = ipi) T , and the condensate wavefunctions 
and ipi are related to the two pseudo-spin components of the 
BEC. The single particle Hamiltonian Ho in Eq. (Q]) reads: 

= ^-(Pxt - k L a z ) 2 + V tr (x)t + Sla x , (2) 
2m 

where p x = —ihd x is the momentum operator in the longitu- 
dinal direction, m is the atomic mass, a XiZ are the usual 2x2 
Pauli matrices, 1 is the unit matrix, kh is the wavenumber 
of the Raman laser which couples the two atomic hyperfine 
states, Vt — \[2VLr is the strength of the Raman coupling, 
while V\, T {x) = muj 2 x 2 /2 is the harmonic trapping poten- 
tial. Finally, the effective ID coupling constants in Eq. (Q}, 
gij = 2aijhuj± = 1, 2), are defined by the s-wave scat- 
tering lengths a.ij \ for attractive interactions, ct^ < 0. 

Let us measure length in units of the transverse harmonic 
oscillator length a±_ = \/h/(muj_\_), energy in units of 
Huj±, and densities in units of 2\au |; furthermore, employing 
the gauge transformation ip^^(x, ^) ~^ ^t,l( x > exp(— z/xt), 
where \i is the chemical potential, we derive from Eq. (Q]) the 
following dimensionless equations of motion for 

idtik = (-^* - ik Ld x + V tt (x) - l^tl 2 - Wll 2 ) tft 
-/i^ t + 0^, (3) 

idt^i = ("^ + ik L d x + Vtr(x) - Wtl 2 - 7l^| 2 ) ^ 
-Hi/jl + Slifo, (4) 

where V tr (x) = (uj x /uj±) 2 x 2 /2, f3 = |ai 2 /an|, 7 = 
|«22/<^n I , and we have used the transformations fc^ — » 
kL/a± and ft — >> QHlj± . 
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Region I: |<l < Q Region II: kt > Q 




FIG. 1 : (Color online) The linear dispersion relation (energy spec- 
trum) uj — uj±(k). The upper branch cj+ has a minimum (/c, uj) = 
(0, in both regions I (left panel) and II (right panel), correspond- 
ing to k\ < Q and k\ > Q. The lower branch cj_ has a minimum 
(maximum) (/c,cj) = (0, — Q) in region I (region II); in region II, 
there also exist two minima (±fco, cj m i n ). 

Limiting cases of the system ©-(|4| with V tr = have been 
studied in a wide range of contexts. First, in the absence of the 
kinetic (ex d 2 ) and self-interaction (|^| 2 ^, IV^I 2 ^) terms, 
the above system becomes the massive Thirring model 111711 , 
which is a Lorentz-invariant completely integrable system of 
classical field theory, possessing exact soliton solutions H 1 811 . 
In the absence of the kinetic terms, but in the presence of self- 
interactions, the same model has been studied in nonlinear op- 
tics; in this case, Eqs. ©-(UJ) take the form of a NLDE, which 
describes solitons in optical fiber gratings A similar 

NLDE was also used in the context of SO-coupled BECs B12I1 
and self-trapped states, in the form of gap solitons, were pro- 
posed. Finally, a model similar to Eqs. ©-Q, which includes 
the kinetic terms with a dispersion coefficient D, was studied 
in Ref. [ 20] ; this model, which finds applications to two cou- 
pled planar nonlinear optical waveguides, supports so-called 
"embedded solitons" for various values of D (and frequency 
uj); these solitons, however, are generally only semi- stable. 

Here, we will use a multiscale expansion method to derive 
approximate soliton solutions of Eqs. ©-(HI) with a frequency 
(chemical potential) residing in the semi-infinite gap of the 
linear spectrum. The solitons will be found to be stable for a 
wide range of experimentally relevant parameter values. Our 
analytical results will be obtained for 7 = 1 and V tr = 0; 
deviations from this choice will be investigated numerically 
and they will not qualitatively alter our results. 

Analytical results. Seeking small- amplitude solutions oc 
ex.p[i(kx — ujt)] of Eqs. ©-© with fi = 0, we obtain the 
following dispersion relation (energy spectrum): 

uj±{k) = ^k 2 ± ^k 2 L k 2 + n 2 , (5) 

which features two distinct branches. The upper branch, 
cj+(fc), always has a minimum at (k,oj) = (0, +Q), and the 
lower branch, u-(k), has different behaviors depending on 
the sign of the parameter A = 1 — k^/Fl: if A > then this 
branch has a minimum (fc,cj) = (0, — Q) (region I), while 
if A < 0, uj-(k) has a maximum (k,u) = (0,—Q) and 
two minima (±&o,cj m m) (region II). The dispersion relation 
uj±(k) is shown in Fig. [TJ clearly, in the linear regime, the 
lowest energy states in region I can only have zero momen- 
tum, k = 0, while in region II they may have either a positive 
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FIG. 2: (Color online) Density profile of the bright soliton in region 
I. Solid line and circles depict, respectively, the analytical result [per- 
taining to Eq. ©] and the numerically found exact solution. Left and 
right insets show, respectively, the real and imaginary parts of the two 
wavef unctions. Parameters are: kh — 8, Q = 120, f3 = 0.8, and 
e 2 uj 0.4. 

or negative momentum, ±fcn, or they can be a linear superpo- 
sition of both modes with momentum ±fco, thus forming the 
"stripe phase" ((J. 

For fi < — ft in region I, or \i < uj m i n in region II, there ex- 
ists a semi-infinite gap where linear modes do not propagate. 
However, matter- wave bright solitons with energies inside the 
semi-infinite gap can be found analytically via a multiscale 
expansion method. In particular, let /i = — Q — c 2 ujo in re- 
gion I and /i = cj m in — e 2 uJo in region II, where e is a for- 
mal small parameter and ujq is a free positive parameter (with 
ujo/Q = O(l)), which sets the energy difference, c 2 ujo, from 
the linear limit inside the semi-infinite gap (cf. Fig. [T]). We 
seek solutions of Eqs. ©-© in the form: 

/ Mx,t)\ _ ( eA(X)\ iKx 

where A(X) and B{X) are unknown functions of the slow 
variable X = ex, while the momentum K is chosen as K = 
in region I and K = ±ko in region II. Expanding A(X) and 
B{X) as a series in e, i.e., A(X) = ^2 n yo^ n a n (X) and 
B{X) = ^Z n>0 e n b n (X), and substituting the above expres- 
sions m Eqs. we obtain the following. 

In Region I, the solvability conditions at the leading [(9(1)] 
and first-order [0(e)] approximations are satisfied if ao = 
-b = u(X) and a x = b x = i(k L /2n)u'(X), where u(X) 
is an unknown complex function (primes denote derivatives 
with respect to X). The latter is determined at the order 
0(e 2 ), where the solvability condition is the following sta- 
tionary nonlinear Schrodinger (NLS) equation: 

u" - \u + v\u\ 2 u = {) 1 (7) 

where the positive coefficients A and v are given by: 

A = 2u; A-\ i/ = 2(l + /3)A- 1 

(recall that A > in region I). 

In region II for K = ±fcn, the solvability condition at the 
leading order reads as a linear equation connecting functions 
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FIG. 3: (Color online) Same as in Fig.[2but for the solitons in region 
II. The top and bottom panels show, respectively, the soliton with k = 
+ko [cf. Eq. J9}] and the "stripe soliton" [cf. Eq. dTTbl . Parameters 
used are as in Fig. but with Q = 35 and e 2 uo = 0.2. 



clq(X) and bo(X), namely 

a = =F k )b = u(X). 

At the first order, we obtain a similar condition for the func- 
tions a\{X) and bi(X), namely 

k L (k L ± h)ai(X) + Qb^X) = i(k L ± k )u'(X). 

Finally, at the order 0(e 2 ), the solvability condition is again a 
stationary NLS of the form of Eq. (|7]), but with the coefficients 
A and v now given by: 



A 



2uj^k\ 
n 



2k L (k L ±k )(ki + k 2 L k$ + f3n 2 ) 



ft 2 /cq 



Taking into regard that the soliton solution of the stationary 
NLS Eq. © is of the form u(X) = yj2\/v sech(VAX), we 
end up with approximate [valid up to the order 0(e 2 )] soli- 
ton solutions of Eqs. ©-dH for the wavefunctions ip^^(x, t). 
These solutions, characterized by the free parameter e^/cJo 
(measuring the energy difference from the linear regime), 
have the following form: in region I, 



and in region II, 



2uj I 2u 

IT^ sech l c V a"* 




f(x) 



-k L (k L ± k ) 



(8) 



,(9) 



where the function f(x) is given by: 

y/2u>okL 



fix) 



sech 




(10) 



Notice that Eq. © describes two different soliton solutions, 
each corresponding to the locations k = ±fco of the energy 
minimum. We can construct still another approximate soli- 
ton solution by using the linear combination of the above ±ko 
soliton states. In particular, Eqs. ©-(HI for 7 = 1 are com- 
patible with the symmetry ^ = —ijj^ (bar denotes complex 
conjugate) and the following solution satisfies this symmetry: 



C\ cos(kox) + iC2 siri(kox) 
—C\ cos(kox) + 1C2 sin(kox) 



(11) 



where C\ = ft + k\ and C2 = — /co&l- It is clear that, op- 
positely to the solutions dU)-® which have a smooth sech 2 - 
shaped density profile, the soliton (ITTb has a spatially mod- 
ulated density profile (with a wavelength 27r//co); thus, this 
"stripe soliton" ([TIT) is directly analogous to the characteristic 
stripe phase of SO-coupled BECs JaLZD, but now for conden- 
sates with attractive interactions. Note that only solutions ([5]) 
and (ITTb were considered in the numerical studies of Ref. [20] ; 
solution © which does not satisfy the symmetry ^ = — ^ 
was not previously explored. 

The above solutions describe different spin polarizations of 
the gas: these are found as the (normalized) longitudinal and 
transverse spin polarization of the solitons, a XyZ = ((T x )/ntou 
where (a XiZ ) = ^a xz ^/ and n to t = l^tl 2 + l^|| 2 is the to- 
tal density. Then, in region I, we find that the solitons are fully 
polarized along the x-axis, i.e., a x = — 1 (and a z = 0). On 
the other hand, in region II, the stripe soliton has again a z = 0, 
while the dzfeo soliton states are characterized by a finite a Z9 
namely a z = =F\/l — (^Ao) 2 an d &x = —ft/k 2 . (with the 
total mean spin being \/cf 2 + cr 2 = 1). Thus, spin polariza- 
tions of the presented solitons bear resemblance to those found 
for nonlinear states in SO-coupled repulsive BECs 0]. 

Stability and Numerical Results. In our numerical simula- 
tions, we have assumed a quasi- ID attractive BEC, confined 
in a harmonic trap with frequencies uj x = 2ir x 20 Hz and 
uj± = 2tt x 1000 Hz containing approximately 10 3 atoms, 
and scattering lengths ratios 1 : 0.8 : 1 (i.e., f3 = 0.8); ad- 
ditionally, we have considered a fixed value of /cl, namely 
&l = 27r/A with A = 804 nm and varied the parameter ft 
in the range (1 -=- 10)E L , with E L = h 2 k 2 L /2m (with m be- 
ing the 7 Li mass), to identify solutions in region I or region 
II (such an investigation complies with pertinent experiments 
with SO-coupled BECs ^\). We have used a fixed-point algo- 
rithm, and an initial ansatz pertaining to solutions (O-© for 
regions I and II, to find respective numerical solutions. Exam- 
ples are provided in Fig. [21 (for region I) and Fig. [3] (for region 
II), where the density profiles, n to t = | 1 2 + 1 I 2 ' as we ^ as 
the real and imaginary parts (insets) are shown; the analytical 
results (solid lines) are in excellent agreement with the nu- 
merical ones (circles and dashed lines). Furthermore, we have 
numerically confirmed (results not shown here) the existence 
of the presented soliton families for 7 7^ 1, in a relatively wide 
range of values, i.e., for 0.5 < 7 < 1.5. 
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FIG. 4: (Color online) Eigenvalues obtained from the spectral prob- 
lem (TT2l) for the stripe soliton (left) and +&o- soliton (right) branches. 
The latter becomes spectrally unstable due to the eigenvalue pair with 
Re(A) 7^ for /i < —42.9. Parameters are as in Fig. [3] 



We have also studied the stability of the solitons. Because 
each solution family corresponds to the energies inside the 
semi-infinite gap, the spectral stability of solitons is controlled 
by the negative index count (explained in Ch. 4 of Ref. 112 ill ). 
Writing the spectral stability problem in the form 




FIG. 5: (Color online) Contour plots showing the evolution of the 
total density for solitons in region I (a) and region II (b)-(d). Panel 
(b) corresponds to stripe soliton, while panels (c) and (d) correspond 
to +/co-solitons with /x = —41.77 > /i c and \i — —43 < /i c re- 
spectively. Other parameters are as in Fig. [2 (but with Q = 70) and 
Fig. [3] but now for 7 = 0.8 and (jo x /(jo z — 0.02. 



Hu = zAJu, (12) 

where u is a 4 x 1 vector of the perturbations to 
$1], H is a 4 x 4 self-adjoint matrix operator 
associated with the right-hand- side of Eqs. ©-© linearized 
around the solitons, J = diag(l, —1,1,-1), and A is a spec- 
tral parameter with the instability growth rate given by Re (A) 
(if positive). The operator H has a finite number of negative 
eigenvalues, denoted by n(H), and a two-dimensional kernel 
spanned by the symmetries of Eqs. ©"©I 

ui = [z^ t ,-^ t ,^,-#J, u 2 = 0o#t>^t>^4>^4J- 

Associated with the eigenvectors of H , there exist generalized 
eigenvectors of the spectral stability problem (Tl2t given by 
solutions of the inhomogeneous equations 

Hvj=iJuj, j = l,2. (13) 

Computing the symmetric matrix of symplectic projections 
with elements = (v^zJiij) = 1,2), where (•,•) 
is a standard inner product, we denote the number of negative 
eigenvalues of D by n(D). The negative index count is now 
given by # = n(H) — n(D) and this number determines the 
number of unstable eigenvalues with Re (A) > and/or the 
number of potentially unstable eigenvalues with Re (A) = 
and negative energy in the spectral stability problem (IT2l) 12111 . 

To assess the stability of our solutions, we have computed 
indices n(H) and n(D) for the soliton solutions in region I 
and II. For solitons in region I and the stripe solitons in region 
II, we have obtained numerically that n(H) = 1 and n(D) = 
1 in their existence intervals, therefore, the negative index # 
is zero. This ensures spectral stability of these solitons. On 
the other hand, for ±/co-solitons in region II, we have obtained 
numerically that n(H) = 3 in the existence interval, but n(D) 
changes from 1 near the bifurcation at /i = u; m i n to 2 for 
smaller values of fi. Therefore, the negative index is # = 2 
near the bifurcation at ji = cj m i n , due to a pair of negative 



energy yet neutrally stable eigenvalues in the spectrum of (fT2l) . 
For smaller values of /i, it switches to # = 1 indicating a real 
unstable eigenvalue. 

These results are confirmed by the numerical approxima- 
tions of eigenvalues in the spectral problem (fT2l) . Figure @] 
shows the eigenvalues associated with the stripe- and +ko- 
solitons in region II. Spectral stability of the former is con- 
trasted with the potential instability of the latter that arises 
when a pair of neutrally stable eigenvalues of negative energy 
crosses zero at /i = /i c « —42.9 and splits along the real axis 
for smaller /i, yielding an exponential growth of perturbations. 

We have also studied the soliton dynamics for 7^1 and in 
the presence of the trap. We have used our fixed point algo- 
rithm to obtain a specific soliton state; then, the numerically 
found soliton was perturbed by a noise of strength « 10% of 
its initial amplitude, and the resulting state was used as ini- 
tial condition for Eqs. O])-© with the parabolic trap. Results 
of direct simulations are shown in Fig. \5\ for 7 = 0.8 and 
trap strength uu x /uu±_ = 0.02. Solitons in region I [panel (a)], 
stripe solitons in region II [panel (b)] and ko -solitons with 
li = —41.77 > /i c , corresponding to their stability region 
[panel (c)], are found to be robust up to t = 4000 (of the order 
of 1 sec in physical units), which was the time of the simula- 
tion. An example of unstable ko -solitons with fi = — 43 < fi c 
is also illustrated [panel (d)] ; in this case, the soliton stays qui- 
escent for small times [see inset in panel (d)], but later starts 
oscillating in the trap due to the onset of the instability. 

We stress that although our analytical results were obtained 
in the case 7 = 1 and Vt Y = 0, the simulations have revealed 
the existence and stability of solitons for a wide range of val- 
ues 7^1, and also in the presence of the trap, as well as for 
different values of f3. This clearly indicates that the presented 
matter-wave soliton families have an excellent chance to be 
observed in experiments with SO-coupled attractive BECs. 

Conclusions. In summary, we have used a multiscale ex- 
pansion method to identify matter- wave bright soliton states 
in SO-coupled BECs with attractive interactions. The soli- 
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tons, which were characterized by a chemical potential resid- 
ing in the semi-infinite gap of the linear spectrum, were found 
in analytical form to exhibit either a smooth (sech 2 -shaped) or 
a modulated density profile, strongly reminiscent of the stripe 
phase of SO-coupled repulsive BECs. Our analytical pre- 
dictions were corroborated by numerical simulations, which 
have shown that the solitons exist and are generally robust for 
a wide range of the physical parameters involved (including 



chemical potential, interatomic interaction strengths and the 
presence of trapping potentials), even in the presence of noise. 
It would be particularly interesting to explore higher dimen- 
sional generalizations of such solitary waves and the potential 
of collapse type phenomenology 1 22] for the various solitonic 
phases discussed above. Naturally, also, experimental realiza- 
tions of such SO-coupled attractive interaction BECs would 
shed considerable light into such investigations. 
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